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Abstract 

A theory is presented for a novel recursion method for 0{N) ab initio tight- 
binding calculations. A long-standing problem of generalizing the recursion 
method to a non-orthogonal basis, which is a crucial step to make the recursion 
method applicable to ab initio calculations, is solved by the introduction of 
the hybrid representation and the two-sided block Lanczos algorithm. As a 
test of efficiency of the new method, convergence properties in energy and 
force of insulators, semiconductors, metals, and molecules are studied for not 
only simple model systems but also some real materials within the density 
functional theory. The present 0{N) method possesses good convergence 
properties for metals as well as insulators. 

The application of the conventional ab initio electronic structure calculations to large 
systems is hampered by their inherent 0(A^ 3 ) scaling properties with N the number of atoms. 
In order to overcome this difficulty, efficient linear scaling algorithms, which are referred to 
as 0(N) methods, have been developed during the last decade |Ip||. A few applications 
of these 0(N) methods have actually demonstrated the power of these 0(N) methods [|J. 
However, several problems still remain in these O(N) methods. In particular, it is well 
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known that the variational O(N) methods produce large errors in the energy of metals ||, 
since in both metals and narrow-gap semiconductors the density matrix p has long range 
correlations in real space compared to that of insulators with a wide gap ||]. Therefore, the 
application of the O(N) methods relying on the locality of p is restricted to materials with 
a wide gap. 

Recently, one of the authors has shown that the block bond-order potential (BOP) is 
a convergent moment-based O(N) method which provides good convergence properties for 
energy and force in metals as well as insulators and semiconductors within an orthogonal 
tight-binding (TB) representation Although it had been well documented that the 
moment-based methods cannot reproduce the vacancy in diamond or silicon within a low 
number of moments 0,0], the block BOP gave the first convergent results for vacancies in 
insulators using a moment-based method with a low number of moments 0. The block 
algorithm guarantees that the a and 7r bonds are treated properly and that the band energy 
is invariant for the rotation of systems ||, and the terminator in the Green functions can 
estimate the long range contributions in a most effective way || for the total energy and force 
calculations. Thus, the block BOP may be applicable to a wide variety of materials with 
reliable accuracy. However, in order to make the method applicable to ab initio calculations, 
we have to remove the limitation of orthogonal basis and reformulate the method in the 
non-orthogonal basis. Several attempts have been made to generalize the recursion method 
to non-orthogonal basis |J. Nevertheless, within our knowledge the formalism is not in a 
satisfactory level and only very few examples of ab initio calculations have been available [JTOj . 



In this Letter we present a new generalized block BOP for the non-orthogonal TB basis by 
introducing a hybrid representation and the two-sided block Lanczos algorithm. We will 
show convergence properties of the method within density functional theory (DFT) and 
demonstrate that the present method is a promising and practical O(N) method in both 
insulators and metals. 

It will be assumed that one-particle wave functions are expanded with a non-orthogonal 
localized basis set {\ia)) where i is a site index, and a an orbital index. Then, the density 



of electrons p(r) in non-spin polarized systems can be written as 

P( r ) = H Xia(r)xjp(r)Q iatj/3 , (1) 

ia,j/3 

where Xia( r ) = ( r Ka), and the bond-order Oiajp is related to the imaginary part of the one 
particle Green function as follows: 

e ia , i/3 = --Im f G iadP (E + i0 + ) f{^f^)dE (2) 

7T J Kb-L 

with the Fermi function f(x) = 1/[1 + exp(x)], where + represents a positive infinitesimal. 
Gia,j/3(Z) is a expectation value of Greenian G(Z) = [Z — H)~ l for dual bases \ia) and \j(5) 
defined by 

\i~a) =J2\j(3)S-l ia , (3) 

where Sjp ia is the inverse of overlap matrix S ia j/3 = (ia\jj3). From Eqs. (1) and (2), we 
see that functionals such as the total energy in DFT can be reformulated as the functionals 
of the bond-order. Therefore, we concentrate on evaluating the bond-order with reasonable 
accuracy and great reduction in the computational efforts. 

Let us introduce a hybrid representation of Hamiltonian which is a non Hermitian matrix 
represented by the original and the dual bases as H' ia -^ = (ia\H\jf3). The hybrid Hamilto- 
nian can be written in the matrix form as H' = S^ 1 !!, where H ia ^p = (ia\H\j/3). With the 
relation G(Z)(ZS — H) = I, the hybrid Green function G'(Z) defined by 

QajfiW = {G(Z)S} ia ^ = (ia\G{Z)\j(3) (4) 

satisfies G'(Z)(Zl — H') = I. One of the merits of using G'(Z) is that its diagonal element 
gives directly the Mulliken population P ia of an orbital \ia): 

Pia = _^i m / G ' iajia (E + + )f(^-f)dE 

= ®ia,j/3SjP,ia- ( 5 ) 
J/3 

In the block BOP, determination of the chemical potential is needed to conserve the total 
number of electrons N e i e in the system , so that the relation of Eq. (5) is very advantageous 
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to computational efficiency because of the simple relation N e i e = J2i a Pia- Thus, we present 
below a prescription how to calculate the hybrid Green functions. The diagonal elements of 
the Green function matrix can be calculated in a numerically stable way by the recursion 
method || based on the Lanczos algorithm [[ll]]. The block BOP method is a general 
recursion method for evaluating efficiently both the diagonal and off-diagonal elements of 
the Green function matrix by the recursion method. Moreover the use of a single site 
containing all the localized orbitals as the starting state in the block Lanczos algorithm 
rather than a single orbital in the usual one conserves the rotational invariance of the total 
energy. In the present case of non-orthogonal basis, we further extend the formalism to 



adopt a two-sided block Lanczos algorithm [12|, since the hybrid Hamiltonian is not any 



more Hermitian. The set of central equations are 

H\U n ) = \U n )A n + \U n ^)B n + \U n+1 )C n+1 , 

{U n \H = A n {U n \ + C n {U n ^\+B n+1 {U n+1 \. (6) 

Am B_ n , and C_ n are recursion block coefficients with x Mj in size, where Mj is the number 
of localized orbitals on the starting atom i, and the underline indicates that the element is 
a block. In the two-sided block Lanczos algorithm the Lanczos vectors in the left and right 
sides have a bi-orthogonality relation. It is essential to start the two-sided block Lanczos 
algorithm with a single site and its corresponding dual state as 

\U ) = (\il),\i2),...,\iM i )), 

|^ ) = (|il),|i2),...,|iAf < )). (7) 

Equation (7) is an optimum choice in terms of computational accuracy and efficiency because 
of the rotational invariance of the total energy and the consistent description for the different 
properties of a, ir, and 6 bonds. 

In the Lanczos basis representation the Hamiltonian H L is block-tridiagonalized as a 
non Hermitian matrix and the Green function matrix G L (Z) is the inverse of the matrix 
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(Zl — H L ), so that the block diagonal element Gq (Z) = (Uq\G\Uq) can be written explicitly 
in the form of the multiple inverse as follows: 

G^ (Z) = [Zl -A - C^Zl -Ai-CA-' ■}' 1 B 2 ]- 1 B 1 }-\ 

(8) 

where the index L indicates the representation based on the Lanczos basis. The off-diagonal 
elements of hybrid Green function matrix can be calculated by using a recurrence relation 
which can be derived basically along the same line as that described for the case of orthogonal 
basis D. The explicit expression consistent with Eqs. (6) and (8) is given below: 



G^ n (Z) = [G^ 1 (Z)(Zl-A n _ 1 ) 

-Gt^{Z)C n _ x - 5 ln l) (B n y\ (9) 



where 5 is Kronecker's delta, and G _ 1 (Z) = (7 = 0. The block elements of the Green 
function matrix have the same relation to the bond-orders based on the Lanczos basis Bq„ 
as Eq. (2) of the dual basis representation. Therefore, we can obtain the bond orders of 
Eq. (2) through the following transformation: 

n,k 

where U_ nj is defined by U_ nj = (U n \(\jl), \ j2), . . . , \jMj}). As a result of the simple inverse 
transformation Eq. (10), we only have to perform the evaluation and the integration of 
the Green functions of the 0th block line in the Lanczos basis representation, which means 



that the computational time of the algorithm is about two times |13| longer compared 
to that of the orthogonal case [0]. Only the hybrid representation can provide this simple 
relation Eq. (10) as well as Eq. (5), while the other representations suffer from computational 
inefficiency PJT0[|. In the generalized block BOP using the non-orthogonal basis we need to 
calculate S 1-1 , the inverse of the overlap matrix. In the following calculations, we used a 



new O(N) efficient method for inverting the overlap matrix [[14 
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In Fig. 1 we show convergence properties of the band energy in an insulator and a metal 
described by a simple s-valent TB as a test of the present method. The errors in the band 
energy at the seven-shell cluster and recursion levels are 0.2 % and 0.9 % for the insulator and 
the metal, respectively. Thus, we see that the block BOP gives sufficient convergent results in 
both the simple insulator and metal. Figures 2(a) and (b) show the error in the band energy 
at the five-shell cluster and recursion levels for insulators and metals described by a simple 
s-valent TB as a function of direct band gap and electronic temperature, respectively. In 
insulators the error goes to zero as the gap increases, while the errors, whose absolute values 
are no more than 0.5 % compared to the band energy in the whole region, are relatively 
small. In metals the error becomes almost negligible for the higher electronic temperature. 
This behavior in both insulators and metals is consistent with the recent study about the 
locality of the density matrix ||, though the block BOP depends on the convergence of the 
moment expansions for the density matrix rather than the locality of the density matrix 
0. From the comparison in the NaCl and FCC structures it is clear that the use of the 
terminator in the diagonal Green functions effectively reduces the error in both cases. 

Next we discuss convergence properties of the block BOP in realistic materials within 



the TB based DFT proposed by Sankey and Niklewski |15||. Figure 3 shows the convergence 



properties of the cohesive energy for carbon in the diamond structure, silicon in the diamond 



structure, fee aluminum, and Cqq molecule [16|. In carbon and silicon the cohesive energies 



rapidly converge to the k-space results in the five and seven-shell clusters, while the conver- 
gence values for the three-shell cluster are in error by 0.4 and 0.9 % from the k-space results, 
respectively. Even for metallic aluminum, the convergence is very fast with respect to the 
number of recursion levels and the errors in the converged values are only 0.3 and 0.1 % for 
the three- and five shell clusters, respectively. For Cqq the convergence is achieved with the 
three-shell cluster. The error at the sixth recursion level is only 0.02 %. As a test of the 
consistency between the total energy and the forces, a constant energy molecular dynamics 
simulations have been performed for diamond within DFT. In Fig. 4 we show the energy 
for diamond at 300 K using five and ten recursion levels. For five and ten recursion levels 
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we see that the total energy is almost conserved, while the ten recursion level calculation 
gives a better result. These test calculations indicate that the block BOP can give forces 
consistent with the total energy. 

In summary, we have presented a theory of the block BOP for O(N) ab initio TB 
calculations. The introduction of the hybrid representation and the two-sided block Lanczos 
algorithm enables us to generalize the theory to a non-orthogonal basis set in a natural way. 
The test calculations for the simple s-valent TB systems and some real systems within DFT 
suggest that the O(iV) method provides a rapid convergence properties for metals as well as 
insulators with sufficient accuracy. Thus, we conclude that the block BOP is a robust O(N) 
methods which is applicable to a wide variety of materials in the ab initio TB approach. 

We would like to thank Y. Morikawa and H. Kino for helpful suggestions about the DFT 
calculations. 
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FIGURES 

FIG. 1. The error, with respect to the standard k-space calculations, in the band energy for an 
insulator (zinc blende) and a metal (FCC) described by a simple s-valent TB model in which the 
nearest neighbor hopping and overlap integrals are -1.0 eV and 0.1, respectively, with others being 
zero, and the number of electrons is the same as that of atoms. The zinc blende has a direct gap 
of 1.0 eV which was controlled by the gap of the on-site energies of the different atoms. In these 
calculations, the seven-shell cluster and a square-root terminator were used. 

FIG. 2. The error in the band energy for (a) insulators and (b) metals, calculated at the 
five-shell cluster and recursion levels. The calculations were carried out with the same s-valent TB 
model as that in Fig. 1 using a square-root terminator. For NaCl and FCC the non terminator 
results are also shown. 

FIG. 3. The error in the cohesive energy for carbon in the diamond structure, silicon in the 
diamond structure, fee aluminum, and Ceo for three-, five-, and seven-shell clusters, calculated 
using a square-root terminator. These calculations were performed within DFT. 

FIG. 4. The potential, kinetic, and total energies as a function of time for molecular dynamics 
simulations of carbon using a three shell cluster and a square root terminator. In panels (a) and 
(b) the results are for five and ten recursion levels at 300 K, respectively. The time step is 0.5 fs. 
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